Datos de estudio

En este ejemplo práctico utilizaremos datos de un estudio sobre la expresión génica en células infectadas con el virus de la gripe (Bajwa et al., 2016). En este estudio se infectó células dendríticas plasmocitoides humanas con el virus de la gripe y se comparó la expresión génica en estas células con la expresión génica en células no infectadas. El objetivo era ver cómo el virus de la gripe afectaba la función de estas células del sistema inmunitario.

Organización de las prácticas

En esta práctica vamos a usar el entorno de desarrollo integrado RStudio (IDE) para R.
R es un lenguaje de programación para computación estadística y gráficos.

El documento actual en el que estamos trabajando es un documento R Markdown. Son totalmente reproducibles y permiten combinar texto, imágenes y código —en lenguaje R.

Para renderizar un documento R Markdown a un archivo HTML, solo necesitas hacer clic en el botón Knit que verás en la barra de RStudio. Este archivo HTML puede compartirse como un informe.

No necesitas renderizar todo el documento cada vez que quieras ver el resultado de tu código en R. Puedes hacer clic en el botón Run current chunk o usar el atajo de teclado Ctrl+Alt+C, y el resultado del código aparecerá debajo.

A lo largo del documento verás diferentes iconos, cuyo significado es:

: información adicional o útil
: un ejemplo
: un ejercicio práctico
: un espacio para responder al ejercicio
: una pista para resolver un ejercicio

Instalación de herramientas

Recomendamos encarecidamente usar un sistema operativo Linux y acostumbrarse a trabajar con la terminal.
En esta práctica, sin embargo, no la usaremos y casi no hay diferencias entre usar RStudio en Windows o en Linux.

Instalación de R

El lenguaje de programación R puede descargarse desde este enlace, y está disponible para Windows, Linux y macOS.

Instalación de RStudio

RStudio se puede descargar desde este enlace.
Se instala fácilmente en Windows, Linux y macOS.

Instalación de un paquete

Los paquetes de R son colecciones de funciones y/o datos desarrollados por la comunidad.

Para instalar un paquete, usamos la función install.packages(), indicando entre comillas el nombre del paquete que queremos instalar.

Instalación de ggplot2

ggplot2 es un paquete de visualización de datos para el lenguaje de programación R (Wickham 2009).
Fue creado por Hadley Wickham, implementando la Grammar of Graphics de Leland Wilkinson —un esquema general para la visualización de datos que divide los gráficos en componentes semánticos (Wilkinson 2010).

Para instalar el paquete ggplot2, utiliza lo siguiente:

install.packages("ggplot2")

Cargar un paquete

Para cargar un paquete, usamos la función library(), indicando el nombre de la librería que queremos cargar.

En el caso de ggplot2, usarás:

library("ggplot2")

¡Ya estás listo para empezar a usar R!

Introducción a R y RStudio

Directorio de trabajo

El directorio de trabajo es la carpeta de tu ordenador donde R busca por defecto los archivos que necesitas importar y donde guardará los resultados que exportes (tablas, gráficos, etc.).

Para saber cuál es el directorio de trabajo (working directory) utilizamos la funcion getwd():

getwd()
## [1] "C:/Users/Marta/OneDrive - UAB/Bioestadística/Prácticas/P1"

Podemos definir una nueva ruta de trabajo con setwd():

setwd("/ruta/nou/directori")

O de manera gráfica dentro de RStudio a Session > Set Working Directory > Choose Directory…).

Actividad 1.1

Decide una ruta de trabajo en tu ordenador con setwd() o de manera gráfica. Comprueba que sea correcta con getwd().

Por ejemplo, puedes escoger la carpeta donde te hayas descargado el material del Moodle.

setwd("C:/Users/eduar/OneDrive - UAB/Bioestadística/Prácticas/P1")
getwd()

En este directorio de trabajo que hayas escogido, realizaremos todos los análisis.

Leer el archivo con los datos de expresión

Lee en R el archivo “GSE68849_expression.csv” que encontrarás en la carpeta data. Puedes abrirlo primero con Excel o con otro programa de visualización de datos (como un editor de texto) para comprobar si tiene nombres de columnas y cómo están separadas.

Actividad 1.2

¿Cuál será el código en R para abrir el archivo? Guarda el contenido del archivo en un objeto llamado expressionData.

Puedes utilizar los comandos que vimos en la diapositiva 12 de teoría del Tema 2: Introducción a R. Ten en cuenta que quizás necesites cargar una librería, como library(readr).

library(readr)

expressionData <- read_csv("data/GSE68849_expression.csv", col_names = T)

Actividad 1.3

Explora el archivo. Utiliza comandos como head(expressionData) y View(expressionData) para explorar el contenido de expressionData. ¿Qué información contiene?

head(expressionData)
## # A tibble: 6 x 4
##   subject    treatment gene  expression_level
##   <chr>      <chr>     <chr>            <dbl>
## 1 GSM1684095 control   IFNA5             83.1
## 2 GSM1684096 influenza IFNA5          10096. 
## 3 GSM1684097 control   IFNA5             97.8
## 4 GSM1684098 influenza IFNA5           8181. 
## 5 GSM1684099 control   IFNA5             81.7
## 6 GSM1684100 influenza IFNA5           7055.
# Contiene los datos de expresión génica en dos grupos de sujetos: control y con infección

Tipos de variables

Actividad 1.4

Para saber el tipo de variable que contiene expressionData, podemos utilizar la función str(). Rellena la siguiente tabla indicando qué tipo de variable contiene expressionData.

Variable Tipo de variable (Cuantitativa / Cualitativa) Tipo específico (Discreta / Continua / Nominal / Ordinal) Observaciones
Subject Cualitativa Nominal Identificador único de cada sujeto
Treatment Cualitativa Nominal Dos categorías: “control” e “influenza”
Gene Cualitativa Nominal Nombre del gen analizado
Expression level Cuantitativa Continua Niveles de expresión génica

Estadística descriptiva

R contiene muchas funciones built-in para realizar estadística descriptiva básica.

Actividad 1.5

Dada la columna con el nivel de expresión de los genes expressionData$expression_level.

Calcula la media, la mediana y la desviación estándar del vector expressionData$expression_level.

mean(expressionData$expression_level)
## [1] 3947.509
median(expressionData$expression_level)
## [1] 131.8218
sd(expressionData$expression_level)
## [1] 7210.434

Tidyverse

Actividad 1.6

¿Cuántos genes distintos hay en el set de datos?

Primero, podemos seleccionar la columna gene de nuestro dataframe expressionData y luego usar distinct() para quedarnos solo con los valores únicos.
Después, aplicamos count() para contar cuántos elementos hay en total.

Como queremos contar todas las filas resultantes, no hace falta especificar ningún argumento dentro de count().

Recuerda usar %>% para concatenar comandos.

expressionData %>%
  select(gene) %>%
  distinct() %>%
  count()
## # A tibble: 1 x 1
##       n
##   <int>
## 1    10

Alternativamente, te mostramos una función que no se ha visto en la teoría y que permite hacer lo mismo: n_distinct junto con summarize:

expressionData %>%
  summarize(n=n_distinct(gene))
## # A tibble: 1 x 1
##       n
##   <int>
## 1    10

Actividad 1.7

Una de las columnas es subject. Extrae todos los datos correspondientes al sujeto GSM1684096 y guárdalos en un nuevo objeto llamado expressionData_GSM1684096:

Utiliza la función filter()

expressionData_GSM1684096 <- expressionData %>% filter(subject=="GSM1684096")

Actividad 1.8

Explora expressionData. Estamos interesados en conocer el nivel de expresión del gen IFNA5 en condición de infección para el sujeto GSM1684096. ¿Cómo lo harías? Combina filter() y select():

expressionData_GSM1684096 %>% filter(gene=="IFNA5", treatment=="influenza") %>% select(expression_level)
## # A tibble: 1 x 1
##   expression_level
##              <dbl>
## 1           10096.

Actividad 1.9

Calcula la media de expresión para cada tratamiento. Deberás combinar group_by() y summarize(). La diapositiva 42 puede servir de guía:

expressionData %>% group_by(treatment) %>%
  summarize(avg = mean(expression_level, na.rm = T))
## # A tibble: 2 x 2
##   treatment   avg
##   <chr>     <dbl>
## 1 control    156.
## 2 influenza 7739.

Actividad 1.10

Crea una nueva columna basada en el nivel de expresión de los genes: <10 → “bajo”,10-100 → “moderado”, 100-1000 → “alto”, >1000 → “muy alto”

En la diapositiva 32 viste cómo hacerlo con ifelse() para dos grupos. Ahora lo haremos con case_when(), que permite trabajar con múltiples condiciones. Un ejemplo que puede servirte:

surv_raw <- surv_raw %>%
  mutate(age_group = case_when(
    age < 12 ~ "infante",
    age >= 12 & age < 18 ~ "adolescente",
    age >= 18 & age < 65 ~ "adulto",
    age >= 65 ~ "anciano"
  ))
expressionData <- expressionData %>%
  mutate(expression_level_group = case_when(
    expression_level < 10 ~ "bajo",
    expression_level >= 10 & expression_level < 100 ~ "moderado",
    expression_level >= 100 & expression_level < 1000 ~ "alto",
    expression_level >= 1000 ~ "muy alto"
  ))

Actividad 1.11

En qué grupo de tratamiento estan los genes con categoría de expresión muy alto? en los controles, en los de influenza, o ambos?

expressionData %>% filter(expression_level_group=="muy alto") %>% distinct(treatment)
## # A tibble: 1 x 1
##   treatment
##   <chr>    
## 1 influenza

Visualización de datos

Actividad 1.12

Con esta actividad aprenderás a hacer paso a paso un heatmap.

Haremos un heatmap (con geom_tile) con el número de genes (gene) por categoría de subject usando los datos del archivo expressionData.

En la siguiente imagen puedes ver el gráfico final que haremos:

Paso 1. Visualizar los datos

Primero, crearemos el objeto del heatmap con la función ggplot, y generaremos un gráfico. Queremos un heatmap donde se muestren los diferentes sujetos en el eje x (subject) y los genes en el eje y (gene), y el color de cada celda refleje cuánto se expresa cada gen dentro de cada sujeto (fill = log(expression_level)). Para ello, utilizaremos geom_tile().

Modifica los ??? del siguiente código:

exp.heatmap <- ggplot(data = expressionData,
                      mapping = aes(x = subject,
                                    y = gene,
                                    fill = log(expression_level))) +
  geom_tile()
exp.heatmap

Aunque es un buen comienzo, podemos mejorar varias cosas:

  • Los ejes se podrían visualizar mejor, ya que no podemos leer los nombres de los sujetos.
  • Estaría bien saber qué muestras corresponden a control y cuáles a infección.
  • Usar otra escala de colores.

Paso 2. Mejorar el gráfico

Para mejorar los ejes, haremos lo siguiente:

  1. Utilizaremos una etiqueta más clara para el título del eje x (por ejemplo, “Sujeto”).
  2. Giraremos las etiquetas del eje x con un ángulo de 45 grados para poder leer los nombres de los sujetos.
  3. Omitiremos el título del eje y.
exp.heatmap <- exp.heatmap +
  labs(x = "Subjecte", y = "") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) # Rotar les etiquetes de l'eix x
exp.heatmap

A continuación, necesitamos separar las muestras de control de las muestras infectadas con el virus de la gripe. Para hacerlo, utilizaremos la función facet_grid() de ggplot, usando la variable treatment para separar los dos grupos.

# Faceting
exp.heatmap <-  exp.heatmap +
  # facet_grid makes two panels, one for control, one for flu:
  facet_grid(. ~ treatment, scales = "free_x", space = "free_x")
exp.heatmap

Añade scales = "free_x", space = "free_x" dentro de facet_grid() para eliminar los espacios vacíos del gráfico. ¡Prueba a no añadirlos y observa qué pasa!

Y finalmente, podemos cambiar la escala de color; por ejemplo, podríamos utilizar la escala de color scale_fill_viridis():

library(viridis)
# Change scale colour
exp.heatmap <- exp.heatmap+ 
  scale_fill_viridis()
exp.heatmap

Si has seguido todos los pasos, ¡ahora tu gráfico debería ser exactamente como el que se proporcionó al inicio!

¿Qué genes están sobreexpresados en el grupo infectado?

  Respuesta:

  IFNW1, IFNA5, IFNA2, IFNA16, IFNA13

Entrega

Entrega en el Campus Virtual el documento renderizado (HTML) (botón knitr). Recuerda eliminar eval = F de los chunks para poder visualizar el resultado del código.

Referencias

Este ejemplo se ha creado basándose en el tutorial de Jeff Oliver para hacer el heatmap.